#include "Adetector.h"
void GetFinalState_detector_jet(CDraw &para, AnalyseClass &analyse, AObject &object, TClonesArray *branchJet, MyPlots *plots, bool &Jjet, bool &Jtagjet){
	bool Jbjet=true, Jcjet=true, Jqjet=true, Jtaujet=true, Juntagjet=true;
	if(para.jet.jets.SWITCH){
		para.debug.Message(4,23,"open jet switch, begin get jet collection");
		GetFinalState_detector_jet_collect(para,analyse, analyse.delphes.branchJet,plots,object.Vjet, Jjet);
		if(!Jjet){
			para.debug.Message(4,23,"not find enough jets, return");
			return;
		}
	}
	else{
		para.debug.Message(4,23,"not open jet switch, get back from jet function");
		Jjet=true;
		Jtagjet=true;
		return;
	}

	if(Jjet){
		para.debug.Message(4,23,"find enough jets, sort jets with their PT");
		for(int i=0;i<object.Vjet.size();i++){
			object.Ljet.push_back(object.Vjet[i]->P4());
		}
		std::sort(object.Ljet.begin(),object.Ljet.end(),TLV_compare_Pt);
		object.Vqjet=object.Vjet;
		object.Lqjet=object.Ljet;
	}


	// Analyse tagjets
	if(Jjet && (para.jet.bjet.SWITCH||para.jet.cjet.SWITCH||para.jet.qjet.SWITCH||para.jet.taujet.SWITCH)){
		para.debug.Message(4,23,"get tagged jet");
		GetFinalState_detector_jet_classify(para,analyse,object.Vjet,object.Vbjet,object.Vcjet,object.Vqjet,object.Vtaujet, object.Vuntagjet, plots, Jbjet,Jcjet,Jtaujet,Jqjet, Juntagjet);
		if(Jbjet&&Jcjet&&Jtaujet&&Jqjet&&Juntagjet){
			para.debug.Message(4,23,"get enough tagged jets");
			Jtagjet=true;
		}
		else{
			Jtagjet=false;
			return;
		}
	}
	else{
		para.debug.Message(4,23,"don't need tagged jets, begin making jet signal");
		Jtagjet=true;
		Jjet=GetFinalState_detector_make_signal_jet(para, object);
		return;
	}

	if(Jbjet){
		para.debug.Message(4,23,"sort bjets");
		for(int i=0;i<object.Vbjet.size();i++){
			object.Lbjet.push_back(object.Vbjet[i]->P4());
		}
		std::sort(object.Lbjet.begin(),object.Lbjet.end(),TLV_compare_Pt);
	}
	if(Jcjet){
		para.debug.Message(4,23,"sort cjets");
		for(int i=0;i<object.Vcjet.size();i++){
			object.Lcjet.push_back(object.Vcjet[i]->P4());
		}
		std::sort(object.Lcjet.begin(),object.Lcjet.end(),TLV_compare_Pt);
	}
	if(Jqjet){
		para.debug.Message(4,23,"sort light jets");
		for(int i=0;i<object.Vqjet.size();i++){
			object.Lqjet.push_back(object.Vqjet[i]->P4());
		}
		std::sort(object.Lqjet.begin(),object.Lqjet.end(),TLV_compare_Pt);
	}
	if(Jtaujet){
		para.debug.Message(4,23,"sort tau jets");
		for(int i=0;i<object.Vtaujet.size();i++){
			object.Ltaujet.push_back(object.Vtaujet[i]->P4());
		}
		std::sort(object.Ltaujet.begin(),object.Ltaujet.end(),TLV_compare_Pt);
	}
	
	para.debug.Message(4,23,"after prepare all tag jets, begin making jet signal");
	Jjet=GetFinalState_detector_make_signal_jet(para, object);
}

void GetFinalState_detector_jet_collect(CDraw &para, AnalyseClass &analyse, TClonesArray *branchJet, MyPlots *plots,std::vector<Jet*> &jet, bool &Jjet){
	int loop1=0,loopj=0,num;
	float pt,y,mass;

	if (para.signal.NUM_jets==0){
		Jjet=true;
		return;
	} 

	para.debug.Message(4,23,"in the detector jet collect: begin");
	num = branchJet->GetEntriesFast();
	AFill_Plot(para,plots->fJetnum,num);

	std::vector<Jet*> jet_tmp;
	for(loop1=0;loop1<num;loop1++){
		Jet *jet_ptr = (Jet*) branchJet->At(loop1);
		pt = jet_ptr->PT;
		int taujet = jet_ptr->TauTag;
		y  = std::abs(jet_ptr->Eta);
		mass = jet_ptr->Mass;
		//if( para.jet.jets.JCUT_pt  [0] == true && ( pt     < para.jet.jets.CUT_ptd  [0] || pt     > para.jet.jets.CUT_ptu  [0] )){para.debug.ShowMessage(5,23,"false: jets single jet pt",  pt  );continue;}//endif
		if( para.jet.jets.JCUT_pt  [0] == true && ( pt     < 20 || pt     > para.jet.jets.CUT_ptu  [0] )){para.debug.Message(5,23,"false: jets single jet pt",  pt  );continue;}//endif
		if( para.jet.jets.JCUT_mass[0] == true && ( mass   < para.jet.jets.CUT_massd[0] || mass   > para.jet.jets.CUT_massu[0] )){para.debug.Message(5,23,"false: jets single jet eta", y   );continue;}//endif
		if( para.jet.jets.JCUT_eta [0] == true && ( y      < para.jet.jets.CUT_etad [0] || y      > para.jet.jets.CUT_etau [0] )){para.debug.Message(5,23,"false: jets single jet mass",mass);continue;}//endif
		loopj++;
		jet_tmp.push_back(jet_ptr);
	}//endfor jet_Cindex

	if(loopj>=1&&loopj>=para.signal.NUM_jets){
		para.debug.Message(4,23,"jet number larger than we need, begin check exclusive cuts.");
		for(loop1=1;loop1<loopj+1;loop1++){
			if(para.jet.jets.CUT_num>1){
				pt   = jet_tmp[loop1-1]->PT;
				y    = std::abs(jet_tmp[loop1-1]->Eta);
				mass = jet_tmp[loop1-1]->Mass;
				if( para.jet.jets.JCUT_pt  [loop1] == true && ( pt     < para.jet.jets.CUT_ptd  [loop1] || pt     > para.jet.jets.CUT_ptu  [loop1] )){
					para.debug.Message(5,23,"false: jets single jet pt  ", loop1,  pt);
					Jjet=false;
					return;
				}//endif
				if( para.jet.jets.JCUT_mass[loop1] == true && ( mass   < para.jet.jets.CUT_massd[loop1] || mass   > para.jet.jets.CUT_massu[loop1] )){
					para.debug.Message(5,23,"false: jets single jet mass", loop1,   y);
					Jjet=false;
					return;
				}//endif
				if( para.jet.jets.JCUT_eta [loop1] == true && ( y      < para.jet.jets.CUT_etad [loop1] || y      > para.jet.jets.CUT_etau [loop1] )){
					para.debug.Message(5,23,"false: jets single jet eta ", loop1,mass);
					Jjet=false;
					return;
				}//endif
			}
			jet.push_back(jet_tmp[loop1-1]);
		}
		if(jet.size()>=para.signal.NUM_jets){
			para.debug.Message(4,23,"find enough jets", jet.size());
			Jjet=true;
			analyse.counter.getCounter("basic_jet");
		}
		else{
			para.debug.Message(4,23,"not find enough jets", jet.size());
			Jjet=false;
			return;
		}

		para.debug.Message(4,23,"fill jet pt and eta into plots.");
		for(loop1=0;loop1<loopj;loop1++){
			if(loop1<para.signal.NUM_jets){
				AFill_Plot(para,plots->fJetPT[loop1],jet[loop1]->PT);
				AFill_Plot(para,plots->fJetY[loop1] ,jet[loop1]->Eta);
			} 
		}
		return;
	}
	else{
		para.debug.Message(4,23,"false: jet num is small",loopj); 
		Jjet=false;
		return;
	}

}





void GetFinalState_detector_jet_classify(CDraw &para, AnalyseClass& analyse, std::vector<Jet*> &jet,std::vector<Jet*> &bjet, std::vector<Jet*> &cjet, std::vector<Jet*> &qjet, std::vector<Jet*>& taujet, std::vector<Jet*>& untagjet, MyPlots *plots, bool &Jbjet,bool &Jcjet,bool &Jtaujet,bool &Jqjet, bool &Juntagjet){
	int loop1=0,loopb=0,loopc=0,looptau=0,loopq=0,loopuntag=0,jnum;
	float pt,mass,y;

	para.debug.Message(4,24,"in the detector tagged jet cut function");
	jnum = jet.size();

	if (para.signal.NUM_jets==0){
		para.debug.Message(4,24,"para.signal.NUM_jets=0 in Gettagjet");
		return;
	}

	if (jnum==0){
		para.debug.Message(4,24,"jet size=0 in Gettagjet");
		return;
	}

	std::vector<Jet*> bjet_tmp, cjet_tmp, taujet_tmp, qjet_tmp, untagjet_tmp;
	para.debug.Message(4,24,"searching for tag-jet candidates.");
	for(loop1=0;loop1<jnum;loop1++){

		pt   = jet.at(loop1)->PT;
		y    = std::abs(jet.at(loop1)->Eta);
		mass = jet.at(loop1)->Mass;

		if(jet.at(loop1)->BTag==1 &&  para.jet.bjet.SWITCH && para.jet.bjet.CUT_num>=1){
			if( para.jet.bjet.JCUT_pt  [0] == true && ( pt     < para.jet.bjet.CUT_ptd  [0] || pt     > para.jet.bjet.CUT_ptu  [0] )){para.debug.Message(5,26,"false: bjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.bjet.JCUT_mass[0] == true && ( mass   < para.jet.bjet.CUT_massd[0] || mass   > para.jet.bjet.CUT_massu[0] )){para.debug.Message(5,26,"false: bjet single jet mass", loop1,mass);continue;}//endif
			if( para.jet.bjet.JCUT_eta [0] == true && ( y      < para.jet.bjet.CUT_etad [0] || y      > para.jet.bjet.CUT_etau [0] )){para.debug.Message(5,26,"false: bjet single jet eta ", loop1,   y);continue;}//endif
			loopb++;
			bjet_tmp.push_back(jet.at(loop1));
		}

		else if(jet.at(loop1)->BTag==2 && para.jet.cjet.SWITCH && para.jet.cjet.CUT_num>=1){
			if( para.jet.cjet.JCUT_pt  [0] == true && ( pt     < para.jet.cjet.CUT_ptd  [0] || pt     > para.jet.cjet.CUT_ptu  [0] )){para.debug.Message(5,26,"false: cjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.cjet.JCUT_mass[0] == true && ( mass   < para.jet.cjet.CUT_massd[0] || mass   > para.jet.cjet.CUT_massu[0] )){para.debug.Message(5,26,"false: cjet single jet mass", loop1,mass);continue;}//endif
			if( para.jet.cjet.JCUT_eta [0] == true && ( y      < para.jet.cjet.CUT_etad [0] || y      > para.jet.cjet.CUT_etau [0] )){para.debug.Message(5,26,"false: cjet single jet eta ", loop1,   y);continue;}//endif
			loopc++;
			cjet_tmp.push_back(jet.at(loop1));
		}
		else if(jet.at(loop1)->TauTag==1 && para.jet.taujet.SWITCH && para.jet.taujet.CUT_num>=1){
			if( para.jet.taujet.JCUT_pt  [0] == true && ( pt          < para.jet.taujet.CUT_ptd  [0] || pt          > para.jet.taujet.CUT_ptu  [0] )){para.debug.Message(5,26,"false: cjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.taujet.JCUT_mass[0] == true && ( mass        < para.jet.taujet.CUT_massd[0] || mass        > para.jet.taujet.CUT_massu[0] )){para.debug.Message(5,26,"false: cjet single jet mass", loop1,mass);continue;}//endif
			if( para.jet.taujet.JCUT_eta [0] == true && ( y           < para.jet.taujet.CUT_etad [0] || y           > para.jet.taujet.CUT_etau [0] )){para.debug.Message(5,26,"false: cjet single jet eta ", loop1,   y);continue;}//endif
			looptau++;
			taujet_tmp.push_back(jet.at(loop1));
		}
		else if(para.jet.qjet.CUT_num>=1){
			if( para.jet.qjet.JCUT_pt  [0] == true && ( pt     < para.jet.qjet.CUT_ptd  [0] || pt     > para.jet.qjet.CUT_ptu  [0] )){para.debug.Message(5,26,"false: qjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.qjet.JCUT_mass[0] == true && ( mass   < para.jet.qjet.CUT_massd[0] || mass   > para.jet.qjet.CUT_massu[0] )){para.debug.Message(5,26,"false: qjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.qjet.JCUT_eta [0] == true && ( y      < para.jet.qjet.CUT_etad [0] || y      > para.jet.qjet.CUT_etau [0] )){para.debug.Message(5,26,"false: qjet single jet pt  ", loop1,  pt);continue;}//endif
			loopq++;
			qjet_tmp.push_back(jet.at(loop1));
		}
		else{
			if( para.jet.untagjet.JCUT_pt  [0] == true && ( pt     < para.jet.untagjet.CUT_ptd  [0] || pt     > para.jet.untagjet.CUT_ptu  [0] )){para.debug.Message(5,26,"false: untagjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.untagjet.JCUT_mass[0] == true && ( mass   < para.jet.untagjet.CUT_massd[0] || mass   > para.jet.untagjet.CUT_massu[0] )){para.debug.Message(5,26,"false: untagjet single jet pt  ", loop1,  pt);continue;}//endif
			if( para.jet.untagjet.JCUT_eta [0] == true && ( y      < para.jet.untagjet.CUT_etad [0] || y      > para.jet.untagjet.CUT_etau [0] )){para.debug.Message(5,26,"false: untagjet single jet pt  ", loop1,  pt);continue;}//endif
			loopuntag++;
			untagjet_tmp.push_back(jet.at(loop1));
		}
	}	



	if(loopb>=para.signal.NUM_bjet && para.jet.bjet.SWITCH){
		para.debug.Message(4,24,"check bjet cuts");
		for(loop1=1;loop1<loopb+1;loop1++){
			if(para.jet.bjet.CUT_num>1){
				pt   = bjet_tmp[loop1-1]->PT;
				y    = std::abs(bjet_tmp[loop1-1]->Eta);
				mass = bjet_tmp[loop1-1]->Mass;
				if( para.jet.bjet.JCUT_pt  [loop1] == true && ( pt     < para.jet.bjet.CUT_ptd  [loop1] || pt     > para.jet.bjet.CUT_ptu  [loop1] )){
					para.debug.Message(5,24,"false: bjet single jet pt  ", loop1,  pt);
					continue;
				}//endif
				if( para.jet.bjet.JCUT_eta [loop1] == true && ( y      < para.jet.bjet.CUT_etad [loop1] || y      > para.jet.bjet.CUT_etau [loop1] )){
					para.debug.Message(5,24,"false: bjet single jet eta ", loop1,   y);
					continue;
				}//endif
				if( para.jet.bjet.JCUT_mass[loop1] == true && ( mass   < para.jet.bjet.CUT_massd[loop1] || mass   > para.jet.bjet.CUT_massu[loop1] )){
					para.debug.Message(5,24,"false: bjet single jet mass", loop1,mass);
					continue;
				}//endif
			}
			else{
				para.debug.Message(4,24,"bjet cut number < 1", para.jet.bjet.CUT_num);
			}
			bjet.push_back(bjet_tmp[loop1-1]);
		}
		if(bjet.size()>=para.signal.NUM_bjet){
			Jbjet=true;
			analyse.counter.getCounter("basic_tagged_bjet");
			para.debug.Message(4,24,"bjet size > signal. NUM_bjet", bjet.size(), para.signal.NUM_bjet);
		}
		else{
			Jbjet=false;
			ShowMessage(4,"bjet fill",loopb, bjet.size(), para.signal.NUM_bjet);
			return;
		}
		for(loop1=0;loop1<bjet.size();loop1++){
			if(loop1<para.signal.NUM_bjet){
				AFill_Plot(para,plots->fBJetPT[loop1],bjet[loop1]->PT);
				AFill_Plot(para,plots->fBJetY[loop1] ,bjet[loop1]->Eta);
			} 
		}
	}
	else{
		para.debug. Message(4,24,"bjet",loopb);
		Jbjet=false;
		return;
	}

	if(para.signal.NUM_bjet!=0){
		AFill_Plot(para,plots->fBJetnum,bjet.size());
	}

	if(loopc>=para.signal.NUM_cjet && para.jet.cjet.SWITCH  ){
		para.debug.Message(4,24,"check cjet cuts");
		for(loop1=1;loop1<loopc+1;loop1++){
			if(para.jet.cjet.CUT_num>1){
				pt   = cjet_tmp[loop1-1]->PT;
				y    = std::abs(cjet_tmp[loop1-1]->Eta);
				mass = cjet_tmp[loop1-1]->Mass;
				if( para.jet.cjet.JCUT_pt  [loop1] == true && ( pt     < para.jet.cjet.CUT_ptd  [loop1] || pt     > para.jet.cjet.CUT_ptu  [loop1] )){
					para.debug.Message(5,24,"false: cjet single jet pt  ", loop1,  pt);
					continue;
				}//endif
				if( para.jet.cjet.JCUT_eta [loop1] == true && ( y      < para.jet.cjet.CUT_etad [loop1] || y      > para.jet.cjet.CUT_etau [loop1] )){
					para.debug.Message(5,24,"false: cjet single jet eta ", loop1,   y);
					continue;
				}//endif
				if( para.jet.cjet.JCUT_mass[loop1] == true && ( mass   < para.jet.cjet.CUT_massd[loop1] || mass   > para.jet.cjet.CUT_massu[loop1] )){
					para.debug.Message(5,24,"false: cjet single jet mass", loop1,mass);
					continue;
				}//endif
			}
			cjet.push_back(cjet_tmp[loop1-1]);
		}
		if(cjet.size()>=para.signal.NUM_cjet){
			Jcjet=true;
			analyse.counter.getCounter("basic_tagged_cjet");
		}
		else{
			Jcjet=false;
			return;
		}
		for(loop1=0;loop1<cjet.size();loop1++){
			if(loop1<para.signal.NUM_cjet){
				AFill_Plot(para,plots->fCJetPT[loop1],cjet[loop1]->PT);
				AFill_Plot(para,plots->fCJetY[loop1] ,cjet[loop1]->Eta);
			} 
		}
	}
	else{
		para.debug.Message(4,24,"cjet",loopc);
		Jcjet=false;
		return;
	}

	if(para.signal.NUM_cjet!=0){
		AFill_Plot(para,plots->fBJetnum,cjet.size());
	}

	if(looptau>=para.signal.NUM_taujet && para.jet.taujet.SWITCH   ){
		para.debug.Message(4,24,"check tau jet cuts");
		for(loop1=1;loop1<looptau+1;loop1++){
			if(para.jet.taujet.CUT_num>1){
				pt   = taujet_tmp[loop1-1]->PT;
				y    = std::abs(taujet_tmp[loop1-1]->Eta);
				mass = taujet_tmp[loop1-1]->Mass;
				if( para.jet.taujet.JCUT_pt  [loop1] == true && ( pt          < para.jet.taujet.CUT_ptd  [loop1] || pt          > para.jet.taujet.CUT_ptu  [loop1] )){
					para.debug.Message(5,24,"false: taujet single jet pt  ", loop1,  pt);
					continue;
				}//endif
				if( para.jet.taujet.JCUT_eta [loop1] == true && ( y           < para.jet.taujet.CUT_etad [loop1] || y           > para.jet.taujet.CUT_etau [loop1] )){
					para.debug.Message(5,24,"false: taujet single jet eta ", loop1,   y);
					continue;
				}//endif
				if( para.jet.taujet.JCUT_mass[loop1] == true && ( mass        < para.jet.taujet.CUT_massd[loop1] || mass        > para.jet.taujet.CUT_massu[loop1] )){
					para.debug.Message(5,24,"false: taujet single jet mass", loop1,mass);
					continue;
				}//endif
			}
			taujet.push_back(taujet_tmp[loop1-1]);
		}
		if(taujet.size()>=para.signal.NUM_taujet){
			Jtaujet=true;
			analyse.counter.getCounter("basic_tagged_taujet");
		}
		else{
			Jtaujet=false;
			return;
		}
		for(loop1=0;loop1<taujet.size();loop1++){
			if(loop1<para.signal.NUM_taujet){
				AFill_Plot(para,plots->fTauJetPT[loop1],taujet[loop1]->PT);
				AFill_Plot(para,plots->fTauJetY[loop1] ,taujet[loop1]->Eta);
			} 
		}
	}
	else{
		para.debug.Message(4,24,"taujet",looptau);
		Jtaujet=false;
		return;
	}

	if(para.signal.NUM_taujet!=0){
		AFill_Plot(para,plots->fTauJetnum,taujet.size());
	}

	if(loopq>=para.signal.NUM_qjet && para.jet.qjet.SWITCH   ){
		para.debug.Message(4,24,"check light jet cuts");
		for(loop1=1;loop1<loopq+1;loop1++){
			if(para.jet.qjet.CUT_num >1){
				pt   = qjet_tmp[loop1-1]->PT;
				y    = std::abs(qjet_tmp[loop1-1]->Eta);
				mass = qjet_tmp[loop1-1]->Mass;
				if( para.jet.qjet.JCUT_pt  [loop1] == true && ( pt     < para.jet.qjet.CUT_ptd  [loop1] || pt     > para.jet.qjet.CUT_ptu  [loop1] )){
					para.debug.Message(5,24,"false: qjet single jet pt  ", loop1,  pt);
					continue;
				}//endif
				if( para.jet.qjet.JCUT_eta [loop1] == true && ( y      < para.jet.qjet.CUT_etad [loop1] || y      > para.jet.qjet.CUT_etau [loop1] )){
					para.debug.Message(5,24,"false: qjet single jet eta ", loop1,   y);
					continue;
				}//endif
				if( para.jet.qjet.JCUT_mass[loop1] == true && ( mass   < para.jet.qjet.CUT_massd[loop1] || mass   > para.jet.qjet.CUT_massu[loop1] )){
					para.debug.Message(5,24,"false: qjet single jet mass", loop1,mass);
					continue;
				}//endif
			}
			qjet.push_back(qjet_tmp[loop1-1]);
		}
		if(qjet.size()>=para.signal.NUM_qjet){
			Jqjet=true;
			analyse.counter.getCounter("basic_tagged_qjet");
		}
		else{
			Jqjet=false;
			return;
		}
		Jqjet=true;
		for(loop1=0;loop1<qjet.size();loop1++){
			if(loop1<para.signal.NUM_qjet){
				AFill_Plot(para,plots->fQJetPT[loop1],qjet[loop1]->PT);
				AFill_Plot(para,plots->fQJetY[loop1] ,qjet[loop1]->Eta);
			} 
		}
	}
	else{
		para.debug.Message(4,24,"qjet",loopq);
		Jqjet=false;
		return;
	}

	if(para.signal.NUM_qjet!=0){
		AFill_Plot(para,plots->fQJetnum,qjet.size());
	}

	return;

}




void GetFinalState_detector_pseudojet(TClonesArray *branchJet, std::vector<fastjet::PseudoJet> &jet){
	int loop1=0,num;
	num = branchJet->GetEntries();


	for(loop1=0;loop1<num;loop1++){
		Jet *jet_ptr = (Jet*) branchJet->At(loop1);
		fastjet::PseudoJet candi=fastjet::PseudoJet(jet_ptr->PT*TMath::Cos(jet_ptr->Phi),jet_ptr->PT*TMath::Sin(jet_ptr->Phi),jet_ptr->PT*TMath::SinH(jet_ptr->Eta),jet_ptr->PT*TMath::CosH(jet_ptr->Eta));
		jet.push_back(candi);
	}//endfor jet_Cindex
}
